# === setup === #
library(here)
source(here("GitControlled/Codes/0_1_ls_packages.R"))
source(here("GitControlled/Codes/0_0_functions.R"))
#/*--------------------------------*/
#' ## Data
#/*--------------------------------*/
# === NRD boundary === #
nrd_boud <-
here("Shared/Data/WaterAnalysis/NRD_bd/BND_NaturalResourceDistricts_DNR.shp") %>%
st_read() %>%
filter(NRD_Name %in% c("Lower Republican", "Tri-Basin")) %>%
st_transform(4269)
Reading layer `BND_NaturalResourceDistricts_DNR' from data source
`/Users/shunkeikakimoto/Dropbox/ResearchProject/HeterogeneousAllocation/Shared/Data/WaterAnalysis/NRD_bd/BND_NaturalResourceDistricts_DNR.shp'
using driver `ESRI Shapefile'
Simple feature collection with 23 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -11583170 ymin: 4865934 xmax: -10609670 ymax: 5312233
Projected CRS: WGS 84 / Pseudo-Mercator
# === well are located in east side of US Highway 183? === #
east_or_west <-
here("Shared/Data/WaterAnalysis/east_or_west.rds") %>%
readRDS()
# === Data period: 2008-2015 === #
reg_data_all <-
here("Shared/Data/WaterAnalysis/comp_reg_dt.rds") %>%
readRDS() %>%
.[,owner_name := paste0(firstname, "_", lastname)] %>%
.[year %in% 2008:2015 & usage <= 45] %>%
east_or_west[., on = "wellid"]
# === Data period: 2008-2012 === #
reg_data <-
reg_data_all %>%
.[year %in% 2008:2012]
# === All Well Data === #
ir_data <-
here('Shared/Data/ir_reg.rds') %>%
readRDS() %>%
data.table() %>%
.[source %in% c('Meter','METER','metered') & nrdname %in% c("Lower Republican", "Tri-Basin")] %>%
.[,owner_name := paste0(firstname, "_", lastname)]
# .[year %in% 2008:2012]
# --- check: Low_Tri_5mi --- #
ggplot() +
geom_sf(data=nrd_boud) +
geom_sf(data=st_as_sf(ir_data, coords = c("longdd","latdd"), crs = 4269),
aes(color=factor(Low_Tri_5mi)), size=0.5)
"Lower Republican ownerA")to filter the data.outside_5mi_LR <-
ir_data %>%
.[,nrd_owner_name := paste0(nrdname, " ", owner_name) ] %>%
.[nrdname=="Lower Republican" & Low_Tri_5mi==0, nrd_owner_name] %>%
unique()
both_InOut_5mi_LR showing whether a landowner has
farmlands in both inside and outside of the buffer facing LR so that we
can distinguish them later.# create "both_InOut_5mi" index
reg_data <-
reg_data %>%
.[,nrd_owner_name := paste0(nrdname, " ", owner_name) ] %>%
.[, both_InOut_5mi_LR := ifelse(nrd_owner_name %in% outside_5mi_LR, 1, 0)]
reg_data %>%
unique(., by="nrd_owner_name") %>%
.[, .N, by=.(nrdname, both_InOut_5mi_LR) ] %>%
.[order(nrdname)] %>%
print()
nrdname both_InOut_5mi_LR N
<char> <num> <int>
1: Lower Republican 0 306
2: Lower Republican 1 81
3: Tri-Basin 0 445
both_InOut_5mi_LR=1 indicates the farmers who also owns
farmland outside of the 5-mile buffer of LR side, and the number fo
those farmers is 81.
both_InOut_5mi_LR to filter the
data.mod_reg_data <-
reg_data[both_InOut_5mi_LR==0,]
# - TB - #
reg_data_TB <-
mod_reg_data[nrdname=="Tri-Basin",]
# - LR - #
reg_data_LR <-
mod_reg_data[nrdname=="Lower Republican",]
tr (combinations of “township” and
“range” index).tr is not one, if a farmer has farmland
across several “range” .# Ideally, a single farmers have only one single "tr". Check how many of farmers whose farmland is scattered across multiple "tr"
tr_cl <-
reg_data_LR %>%
distinct(nrd_owner_name, tr) %>%
.[, index := 1,] %>%
dcast(.,nrd_owner_name ~ tr, value.var = "index") %>%
replace(is.na(.), 0) %>%
.[, total := rowSums(.SD), by = nrd_owner_name]
# tr_cl[total>1, ]
tr
segments.county_cl <-
reg_data_LR %>%
distinct(nrd_owner_name, county) %>%
.[, index := 1,] %>%
dcast(.,nrd_owner_name ~ county, value.var = "index") %>%
replace(is.na(.), 0) %>%
.[, total := rowSums(.SD), by = nrd_owner_name]
# county_cl[total>1, nrd_owner_name]
county_cl[total>1, ] %>% print()
Key: <nrd_owner_name>
nrd_owner_name FRANKLIN FURNAS HARLAN KEARNEY total
<char> <num> <num> <num> <num> <num>
1: Lower Republican David_Black 1 0 1 0 2
2: Lower Republican John & Ingrid_Tangeman 1 1 0 0 2
3: Lower Republican Wayne & Linda_Brummer 1 0 1 0 2
4: Lower Republican _Franklin County Farm 1 0 0 1 2
ne_county_sf <- tigris::counties("Nebraska", cb = TRUE, resolution="20m", progress_bar = FALSE) %>%
dplyr::select(COUNTYFP, NAME) %>%
filter(NAME %in% str_to_sentence(unique(mod_reg_data$county)))%>%
st_transform(4269)
county_cl_sf <-
reg_data_LR %>%
.[nrd_owner_name %in% county_cl[total>1, nrd_owner_name],] %>%
distinct(nrd_owner_name, longdd, latdd) %>%
st_as_sf(., coords = c("longdd","latdd"), crs = 4269)
ggplot()+
geom_sf(data=ne_county_sf, aes(fill=NAME), alpha=0.6) +
geom_sf(data=nrd_boud, color="blue", fill=NA) +
geom_sf(data=county_cl_sf, size=1)+
facet_wrap(~nrd_owner_name, ncol=2)
# reg_data_LR %>%
# .[nrd_owner_name == "Lower Republican John & Ingrid_Tangeman"] %>%
# distinct(nrd_owner_name, longdd, latdd, .keep_all = TRUE)
# reg_data_LR %>%
# .[nrd_owner_name == "Lower Republican David_Black"] %>%
# distinct(nrd_owner_name, longdd, latdd, .keep_all = TRUE)
Each panel in the figure above shows how farmland of each of those four farmers is distributed.
# reg_data_LR$county %>% unique()
reg_data_LR2 <-
reg_data_LR %>%
.[,county_fix := case_when(
# - first group - #
nrd_owner_name == "Lower Republican John & Ingrid_Tangeman" ~ "FURNAS",
nrd_owner_name == "Lower Republican David_Black" ~ "HARLAN",
# - second group - #
nrd_owner_name == "Lower Republican _Franklin County Farm" ~ "FRANKLIN",
nrd_owner_name == "Lower Republican Wayne & Linda_Brummer" ~ "HARLAN",
TRUE ~ county
)]
# county_cl2 <-
# reg_data_LR2 %>%
# distinct(nrd_owner_name, county_fix) %>%
# .[, index := 1,] %>%
# dcast(.,nrd_owner_name ~ county_fix, value.var = "index") %>%
# replace(is.na(.), 0) %>%
# .[, total := rowSums(.SD), by = nrd_owner_name]
# county_cl2[total>1, ]
reg_data_LR_mean_final <-
reg_data_LR2 %>%
.[,.(
usage = mean(usage),
treat2 = mean(treat2),
# --- soil --- #
silttotal_r = mean(silttotal_r),
claytotal_r = mean(claytotal_r),
slope_r = mean(slope_r),
ksat_r = mean(ksat_r),
awc_r = mean(awc_r),
# --- weather --- #
pr_in = mean(pr_in),
pet_in = mean(pet_in),
gdd_in = mean(gdd_in),
# --- tr --- #
county = unique(county_fix)
), by = .(nrd_owner_name, year)] %>%
.[,county_year := paste0(county, "_", year)]
reg_data_TB_final <-
reg_data_TB %>%
.[, names(reg_data_LR_mean_final) %>% .[-length(.)], with=FALSE] %>%
.[,county_year := paste0(county, "_", year)]
reg_data_final <- rbind(reg_data_LR_mean_final, reg_data_TB_final)
county_yearusge is the average of water use (acre-inch).usge by owner
+ So I needed to calculate with weighted average. (weighted by
acres). + We have two i.acres and
acres + Need to have correct data. Can I reproduce the data
from scratch?